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We investigate dynamic properties of inhomogeneous nano-materials, which appear in analytical 
descriptions typically as series of 5-functions with corresponding Gibbs weights. We focus on ob- 
servables relevant for transport theories of Josephson junction arrays and granular systems near the 
superconductor - insulator transition. Furthermore, our description applies to the theory of tunnel 
junctions exchanging energy with a "bath" , the latter having a discrete spectrum. Using the matrix 
theta-function formalism we find an analytical expression for the transport characteristics capturing 
the complete temperature driven transition from the quantum to the classical regime. 
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I. INTRODUCTION 



The investigation of observables having singular or 
highly oscillating behaviors at microscopic scales, which 
become smooth after appropriate coarse graining is a typ- 
ical occurrence in many research areas of modern physics. 
A detailed theoretical investigation and description is 
therefore highly desirable and important. Quasiclassi- 
cal Green's functions have been introduced in the trans- 
port theory of spatially nonuniform superconductors and 
superconducting hybrid structures as an "envelop" ap- 
proximation for exact Green's functions quickly oscillat- 
ing at Fermi wavelength scales, see Ref. [lj for a review. 
One result of this coarse-gaining of the order parameter 
is the Landau-Khalatnikov theory of second order phase 
transitions 0. Envelop approximations of highly oscil- 
lating solutions of hydrodynamic equations gave an op- 
portunity to achieve analytical and numerical advances 
in hydrodynamics Q . In the study of granular supercon- 
ducting materials, Josephson junction arrays (JJa) and 
strongly disordered superconducting films near the su- 
perconductor - insulator transition, the current-voltage 
characteristics and the dynamic conductivity show sin- 
gular behaviors described by weighted superpositions of 
delta-functions (4T-flfj|. The imaginary part of the two- 
particle Green's function in ultrasmall metallic granules, 
where the electron spectrum is essentially discrete, have 
the form of a series of delta-functions. The transport 
characteristics of ultrasmall tunnel junctions exchanging 
energy with a quantum bath with discrete spectrum have 
a similar singular behavior [TT|-il3| . An important and in 
general unsolved problem is the construction of the en- 
velop approximation for an observable (which is typically 
measured) in case the theoretically derived expression is a 
countable superposition of delta-functions. Remarkably, 
this problem has a general analytical solution going be- 
yond standard approximations that allows transforming 
the discrete series into integrals. We found the envelop 
approximation for singular observables analytically, using 
the Jacoby theta-function formalisms flij . Our calcula- 



tions can be physically interpreted using the language of 
the Landau-Hopf turbulence [l(| EH, EH ; a * least quali- 
tatively. 

In this paper we derive the envelop approximation for 
the real-valued function I(u>) with real argument lj: 



/M = -]Tlm 



L P 

1 n 



cu — cu n + i0 



= ^p n ^-^ n ). (i) 



Here, n, is a ./V-dimcnsional vector of quantum numbers, 
w n is a scalar real function of n and the summation is per- 
formed over all configurations of n with weights (proba- 
bilities) P n > 0. In physical applications lj is usually a 
variable with dimension of energy, w n is closely related 
to (the difference of) the energy levels of a quantum sys- 
tem, while P n is the Gibbs distribution probability. In 
the most general situation when all components of n be- 
long to a countable set, standard approximations like 
the Euler-Maclaurin asymptotic [l~7| do not work any- 
more. In that case one can not naively integrate out 
the 5- functions in ([T]). Our main result deals with this 
case. We find the conditions when the density of the 5- 
functions starts to increase similar to the development of 
Landau-Hopf chaos and find the envelope approximation 
for I(lj) in that case. We discuss the use of our results for 
solutions of specific physical problems mentioned above. 

The case with linear w n and P n with bilinear exponent 
is one of the most relevant for applications. We focus on 
the situation when 



LJ n = e • n, 



P n = CXp - 



n T • E ■ n 



(2) 



Here n = (711,712, . ..,Hjv) t is a vector of integer num- 
bers, e is a (fixed) vector and E is a symmetric positively 
defined TV x N matrix. We will also show that our inves- 
tigation is not crucially dependent on the specific form 
of P n in ©. 

The case where P n decays quickly with |n| (all eigen- 
values A; of the matrix E are much larger than unity) 
is the "quantum" limit. Then the set of numbers n is 
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essentially discrete. Moreover P n effectively restricts n 
to a subset with small |n|. In this case the graph of I(ui) 
looks like a sparse sequence of isolated (5-function peaks. 

Rather straight-forward is the opposite "classical" case 
where P n depends only slightly on n (all eigenvalues Ai 
of the matrix E are much smaller than unity and N is 
sufficiently large). Then the Euler-Maclaurin approxima- 
tion is applicable and one can replace the discrete sum 
over n in Eq. (TTJ by an integration over dn. The result 
is J(w) cx exp(— Lu 2 p/2), where p is constant. 

The most interesting case is the "mesoscopic" case 
when the eigenvalues A.; of the matrix E are slightly 
smaller or of the order of unity and iV is not too large. 
This case is in between the quantum and the classical 
case. Then the restriction for the choice of n is rather 
weak and there are a many n- vectors solving the inequal- 
ity 



\ui - e • n < a, 



(3) 



where a <C 1 is the width of the 5-function. [(5-functions 
in physical applications have always some small width, a, 
due to, e.g., interaction with a heat bath (dissipation).] 
If the components of the vector e are integer numbers 
(commensurate) then only if uj is close to an integer [the 
fractional part of uj is smaller than a\ there is vector n 
- the solution of inequality (PJ. However, in case the 
components of vector e are not commensurate - as it is 
most natural in physical systems - the formal solution of 
Eq. ([3]) for arbitrary u> is only applicable in the classical 
regime where no (integer) restriction of |n| exists. In the 
mesoscopic regime, n T En < 1, the set of cj for which 
Eq. ([3]) has a solution is restricted. But the effective mea- 
sure of this set is much larger then in the case when e has 
commensurable components. This implies that the effec- 
tive density of (5-functions in Eq. (ffj strongly increases 
in a given lu interval, but is still discrete. We will fur- 
ther refer to this observation as the <5-function "conden- 
sation" . As the result the singular part of becomes 
relatively small and I(u) has a smooth envelop approx- 
imation, cx exp(— u! 2 p/2) , d(uj). We emphasize that the 
functional behavior of I(uj) is different from the classi- 
cal limit due to the nontrivial factor i?(w). We find that 
i? is the one-dimensional Jacoby theta- function, i?(w,r*), 
where the real parameter t* > depends on E and e and 
qualitatively shows to what degree the components of e 
are incommensurate. We explicitly calculate r*, which 
behaves in the classical limit as r* — > oo resulting in 
-> 1. 

I 



In order to investigate / in the regime where the 6- 
functions condensate, it is quite ineffective to use the 
(5-function representation of /, Eq. (JT]), directly. Instead, 
we rewrite I in terms of matrix 6-function, 6(z, T) [gen- 
eralizing the approach suggested in Ref. [lOj]. where the 
real symmetric matrix T has one zero eigenvalue, while 
the other eigenvalues are positive. This zero mode is the 
manifestation of 5- functions in Eq. fTJ). We investigate 
the class of the O-functions with the (nearly) degenerate 
T; it is an interesting problem itself. Using the ma- 
trix 0-function representation of / we can analytically 
and numerically investigate I(oS) in the regime where 
the (5-function condensation takes place. In a limiting 
case we reproduce the results of Ref. [TJ3]. We relate 
the nature of the strong increase of the (5-function den- 
sity to the chaotic behavior of quasi-periodic functions. 
Our results can help understanding the transport the- 
ory of Joscphson junction arrays and the superconductor- 
insulator transition. Finally, we discuss how stable our 
results related to the (5-function condensation in Eq. (TTj 
are with respect to the choice of the shape of the weight 
functions P n other than given in Eq. ((T|). 

The structure of our paper is the following: In Sec. |TT] 
we show how problem (JTJ) can be reformulated in terms 
of matrix theta-functions; In Sec. IIIII we investigate the 
properties of the matrix thcta-function in the mesoscopic 
regime and in particular formulate the conditions for the 
delta-function condensation; In Appendix [C] we give a 
numerical receipt for an efficient calculation of the matrix 
theta-function; finally in the discussion section IIV1 we 
show that the problem we address in this paper has a 
number of important physical applications. 



II. GENERALIZED MATRIX 
THETA-FUNCTION 



We rewrite the sum in Eqs. fl])-© using the Poisson's 
formula for summation 

E /(»)= E / f(x> 2mmX dx, (4) 

, ^ „ J — oo 



n— — oo 



m— — oo 



where / is a continuous intcgrable function. Function / 
in Eq. ^ depends on the vector n = (m, n 2 , ■ ■ ■ , n^y . 
Therefore, we have to introduce also m = (mi, . . . , m^Y 
and x accordingly. Approximating the (5(z)-function by a 
smeared Gaussian function e~ z / 2<T /\/2it<t 2 with a — > 
we get after the integration over x: 



= E 



v /2 7 rcr 2 (27r) Ar detG 



exp 



2a 2 



. .n 2ituj 
+ — (2tt) rriirrij + i — —[eirrij + ejUii] + 



etej] \{xiXj) \. (5) 



Here 



Gij = Ey + e;e,-/7, 7 = a 1 
and (xiXj) = [G~ ]»_ 



(0) 
(7) 



A. cr-expansion 

We simplify Eq. ([5]) by expanding it over a. Then we 
find, 
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FIG. 1. (Color online) The graph of 0(z,t) 



where the explicit 
lirricr-M) 7 det G, p = 



analytical 
- lim^o 7 



expressions for 

-1U-1„T . Q-l . 



[7 X e' 



fJ = 
e - 11 



and a = lim <J _ ! .o 7 1 G 1 e can be found in Appendix, 
Eqs. (1aT1)-(IA71). 

Introducing the Riemann theta- function [3 EH j 



e( w ,r) = ]Te* 



27rcjm-a— 7rm T 7~m 



we rewrite / as follows: 

cxp(— uj 2 p/2) 



1(C): 

T 



2nK. 
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e(w,r), 



(9) 

(10) 
(11) 



The proof that Te = is given in Appendix, Eq. (|A6|) . 
Therefore, the 0-function matrix parameter is degen- 
erate resulting in the appearance of the (5-functions in 
Eq. flj. It is worth noting that a is generally not an 
eigenvector of T (such that Ta ^ 0). 

An important role play the eigenvalues Ti of matrix 
T. We will distinguish the cases when r, arc smaller or 
larger than unity below. Without loss of the generality 
we can assume that the eigenvalue corresponding to the 
eigenvector e has the index one: t\ = 0. Then for the 
other eigenvalues follows: Ti>i > 0. 

We kept a ~ above which resulted in t\ = 0. In 
typical physical realizations, the <5-function always have 
a finite width a as was mentioned in the introduction. 
If we preserve the sub-leading terms in a > in / then 
Ti ~ a 2 <C Tj>i. Perturbations of the other parameters 
in Eq. Q produce sub-leading corrections to the shape 
of I(cj). 



B. The "minimal model". 

There are special cases when a = e and the parame- 
ters p, g, and T have very simple forms. One case corre- 
sponds to a diagonal matrix E as follows from the proof 
given in Appendix [B] The case, when E has the addi- 
tional structure: Eij = /3e.;<5.y, where j3 is a parameter 
and /3ej > 0, is important for transport physics in JJAs 
and dirty superconducting films near the superconductor- 
insulator transition [f| [H| • We will refer to this case as 



to the "minimal model". It should be noted that e in 
this case is not an eigenvector of E in general. This case 
is of theoretical interest since it is possible to find G^ 1 
exactly and det G analytically in all orders of 7: 



1 



Then we get provided J2i e « 7^ 0: 



1 



= 2n (Sy 
Similarly we find using Eq. (|A3[) that 



P 



P 



(12) 
(13) 

(14) 
(15) 



Expressions Eqs. (fT5|) - (fT5|l generalize similar results in 
Ref. [HI obtained for the limit where j3 > 0, 7 — > 0, 
J2k e k = A, and e is close to the vector (1,1,..., l) r . 

We took a = in Eq. ([14]) such that t\ = 0. If we 
take into account a <^ 1 then we should correct all ma- 
trix elements, Tij, by an additional term, 2ir c 2 /[(X]i e i) 2 ] 
corresponding to t\ = 2ira 2 /\e\ 2 . 



III. MATRIX 6-FUNCTION IN THE 
MESOSCOPIC REGIME 

Below we study how the properties of the function in 
Eq. (fTU|) depend on the parameters e, E, and especially 
N. 



A. N=l 

We start our analysis with the simplest case, N = 1. 
Then we get from Eq. ((HJ) that reduces to the usual 
Jacobi thcta-function: 



0(w,r)=X\ 



,'1'2'KUJTYl — "KTTfl 



(16) 



For t > 1, #(w,t) 1 and for r < 1, #(a;,Ti) sa 
En'K -' — n), see Fig. [1] for an illustration. The case 
= 1 is marginal to some extent because then r = t\ ~ 
a 2 and so 1? is always a set of (5-functions; formally there 
is no classical limit for N = 1. In the following we focus 
on the case N > 1 only. 



B. N > 1 

The classical limit as follows from Eq. (0| can be found 
by setting m = in Eqs. ©-([TO]). Then 0(w,T) -> 1 
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FIG. 2. (Color online) The dots schematically show the set 
m for N = 2. Without loss of the generality we assume that 
e (direction singled out by the green line) is close to e' ' = 
(1, 1). The dotted lines are parallel to e' '. The red points 
belong to the subset of m that gives the leading contribution 
to G in the mesoscopic regime, Eq. @. 



and as follows from Eq. (JTH 



exp(— co 2 p/2) 



N+l, 



(17) 



We recall that in the mesoscopic regime A, <C 1 (toi ^S> 
1). In the mesoscopic (quantum) regime J(w) strongly 
differs from Eq. (fl"7|) because the prefactor, 0(o>,T), in 
Eq. (|10j) behaves in a non-trivial way. 

The definition of the matrix 0-function, Eq. ([9]), in- 
cludes the sum over m. The components of m are inte- 
ger numbers. So m can be treated as edge vectors of the 
nodes of an effective cubic crystal in a iV-dimcnsional 
space, see Fig. [2j It follows from the definition of the 
0-function, Eq. ((9]), that the leading contribution to 
gives m nearly collinear to e. 

We focus first on the case when e is close to the crys- 
tallographic axis e(°) of the effective crystal. Then the 
subset of m, m' ) = me^°\ where m = 0, ±1, . . ., gives 
the leading contribution to 0: 



0(w T) « e i27ro;me ( " ) -a-7rm 2 (e (0) ) T re (0) 
m 

^e(°)-a,T*), (18) 

and 

t* = (e(°') T re (0) = 6e T T6e, (19) 
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FIG. 3. (Color online) The graph of Q(u, T) for N = 8 [black 
color]. The minimal model has been used. Taking the vector 
from the set of the gaussian random numbers with zero aver- 
age and the dispersion a we generated Se = e — e'°\ where 



,(0) 



1, 



1) T . The red line corresponds to $(w,t*) - the 



envelop function (|2ip . The inset shows that the 8- function 
width a is the minimal characteristic frequency of the matrix 
0-function. 



where 5e = e — e'°\ The vectors e and e^ ' should be 
close such that the following condition is fulfilled: 



<n > l, 



(20) 



which ensures a complete delta-function overlap (we re- 
call that the S- function width is oc <r, while t± oc a 2 ). 
Here the integer m max = Int[l/|e — e(°)|], where Int is 
the integer part. This condition indicates how far from 
the origin the distance between the points on the lines 
Oe and Oe^ ^ becomes of the order of the effective crys- 
tal period, see Fig. [5] 

Expression Q18|) is i n f ac t the envelop approxima- 
tion for the matrix theta-function. The ©-function 
quickly oscillates as the function of w on the scale of 
the order of a on top of the envelop function (fT5|) . see 
Fig. [3J The amplitude of the oscillations does not ex- 
ceed exp(— 7TTi(m m ax) 2 )- It is worth noting that con- 
dition (|20p is sufficient but not necessary. Even beyond 
the limitations set by condition (|20l) the envelop function 
Eq. (fTS)) usually still approximates quite well. This 
case is illustrated in Fig. [3] Numerical calculations [l(| 
show that the envelope of the 0-function follows well 
Eq. (JISJ) while |e - e<°)| < ae ANlnN , where A is a con- 
stant of order unity. This exponential factor is closely 
related to the density of the 5- functions in Eq. |T]), which 
have in the mesoscopic regime weights of the same order. 

If T * i5 0-2 then ^(uie^ ■ a, r*) reduces to a super- 
position of (^-functions, see Eq. (fl6|) and Fig. [T] These 
small values of r* appear when e = e'"', see Eq. (|19p . In 
case the components of e are integer numbers (they are 
commensurable) then the (5-function singularities of I(u>) 
remain in the mesoscopic regime as well. However, if e 
and e*- ) are even slightly different, t* can easily become 
of order unity and becomes a smooth function of u). 
Returning to Eq. ((T|) we will have a strong increase of 
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FIG. 4. (Color online) Density graph of 6(w, T) for N = 8. 
The axes OX corresponds to oj while OY axes corresponds to 
P^ 1 . When /3 > 1 the theta function behaves as the discrete 
set of 5-peaks. With decreasing j3, more peaks become visible, 
for j3 ~ 1 the peaks start merging and for /3 <§C 1 the matrix 
theta function becomes nearly smooth. 



the (5-function density and their final overlap, when the 
components of e become not commensurable. 

To conclude this section, we emphasize that Eq. (fT8|) 
is one of our main results. It shows that in the meso- 
scopic regime I(u>) is nearly continuous (plus a relatively 
small and quickly oscillating background) and the form 
of the function is nontrivial. We find the envelop 
approximation of analytically: 



I{uj) - t?(we 



(0) 



exp(— ui 2 p/2) 



N+l 



9 



(21) 



The accuracy of the prefactor in Eq. (|21j) is illustrated 
in Fig. [3] The density graph 2] shows the evolution of 
the matrix O-function in the minimal model when the 
parameter /? switches the model from the quantum to 
the mesoscopic regime. The numerical calculations are 
briefly described in Appendix [C] 



IV. DISCUSSION 

A. Time representation 

The time representation helps to understand the prop- 
erties of I(lu) from another point of view. For a spe- 
cial case it was shown in that the Landau-Hopf tur- 
bulence scenario can account for the merging of the 5- 
f unctions in Eq. ([T]) in the mesoscopic regime. Here we 




FIG. 5. A graphical representation of Eqs. (|22[1 - (|23[ ) describ- 
ing an open winding path on a two-dimensional torus surface. 



apply the ideas developed in Ref. [T(| for the general case. 
To this end we express I (to), Eq. ([1]), in time representa- 
tion, I(t) = J I(uj)e- lut duj, and find: 



I(t)=5>e <EiL * 

11 

tpi(t)=Uit, i 



(22) 
(23) 



Function I(t) belongs to the set of quasi-periodic func- 
tions !19f. If we take ipi as coordinates in an TV- 
dimensional space then the trajectory <fi(t) is modeled 
by a curve on a torus T that wraps around without ever 
exactly coming back on itself if the u>i are incommensu- 
rate, sec Fig. [5] for N = 2. The path covers the torus 
surface densely everywhere. If we return to Eq. ([TJ then 
this property would mean that the <5-function positions 
are densely distributed. In the quantum regime only n 
with |n| < 1 contribute to I(t) such that not more than 
one frequency appears in the exponent in Eq. (|22j) . The 
classical limit formally corresponds to a torus with an 
infinite number of dimensions. 

Our consideration for the condensation of the 8- 
functions in Eq. (JT|) can be extended to the more general 
case of arbitrary weights P n > which decay quickly 
with n [2(|. In this case, the time representation of / 
would consist of quasi-periodic functions as well and the 
topological argument of the path covering a torus surface 
densely would be applicable again. In this more general 
case I (to) cannot be reduced to the matrix theta-function. 
However, the Poisson transformation of / with an appro- 
priate cr-expansion still allows to reduce / to an "easy to 
handle" form for analytical and numerical investigation 
in the regime when the (5-functions start to merge. 

It is worth noting that expressions similar to Eq. (|22[) 
appear in many applications. For example, it describes 
the velocity field v(r) of a liquid when turbulence devel- 
ops according to the Landau-Hopf scenario [l(| EE EH , 

v(t) = J2 A piP2-pn ex P {iJliLiPi </>*(*)}■ Here tnc de- 
veloped turbulence corresponds to the case of large N 
[2l| . Quasi-periodic functions describe quasi-periodic 
motion of mechanical systems [22j and often appear in 
the theory of differential equations, see e.g., Ref. [23j . 
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B. Superconductor-insulator transition in granular 
systems and JJAs 

In the study of granular superconducting materials and 
Josephson junction arrays (JJA) transport properties are 
of central importance. The most interesting case corre- 
sponds to the quantum limit when the Josephson cou- 
pling Ej between the granules is much smaller than the 
characteristic Coulomb energy E c of the Cooper pair lo- 
cated on a granule. The transport problem was solved 
within linear response approaches, see e.g., 0-0- It was 
shown that the conductivity a JJA of the granular su- 
perconductor is a singular function in the leading order 
over Ej. It behaves as the countable superposition of 
(5-functions with the Gibbs weights. 

For example, the conductivity in a JJA according to 
Ref. Q has the following form [which can also be derived 
from Eq.([T]) in appropriate limiting cases]: 

OjM ~ £ - b ini ) exp (- (e * )2 gf tf " 3 ' ) , (24) 

n 

where u> is the frequency, T is the temperature, the in- 
dices i, j label the islands of the JJA, B is the inverse 
capacitance matrix of the JJA and bi = e*(Bu — I?2i)- 
Here e* = 2e is the charge of the Cooper pair and e 
is the electron charge [we will use the system of units 
where e ~ ks = fi = c = 1]. The integer numbers rii 
show the effective number of Cooper pairs sitting on the 
island i. We assume here and below the Einstein sum 
convention. The exponential weights in Eq. (|24|) come 
from the Gibbs distribution. The delta-functions ensure 
energy conservation during the processes of Cooper pair 
tunneling from one granule of the JJA to its neighbor. 
They cannot be easily integrated out in Eq. (|24j) since 
changes of the Coulomb energy are essentially discrete 
because of charge quantization within granules. 

The size of the JJA arrays can be quite small, espe- 
cially in the case of one-dimensional (ID) arrays. There- 
fore, the number of junctions N in the JJA is typically 
finite, N > 10. On the other hand the interaction matrix 
Bij can be rather quickly decaying with \i— j\. Generally, 
we cannot use standard statistical physics approaches 
based on the thermodynamic limit N — > oo trying to 
smear out the <5-function singularities in the observables, 
as it was done in Eq. (|24l) . An important question is 
to understand the nature of the ^-function singularities 
in the observables and finding systematic ways of their 
regularization. 

One way to overcome the difficulties related to the 6- 
functions in transport observables of the JJA was pro- 
posed in Ref. Q. That calculation was based on the as- 
sumption of an energy band for Cooper-pair tunneling. 
However, this band can form in JJa with nearly identi- 
cal granules but should be suppressed by the disorder in 
typical disordered JJAs [l(| ■ 

In Refs. 0, 0, [l(| an attempt was made to find the 
transport characteristics in the disordered JJAs. The 



model of Refs. 0, @ leads to the conductivity behav- 
ing according to Eq. ([T]) with Eij cx E c ei5ij and e close 
to the vector (1,1,..., l) r (a limiting case of our mini- 
mal model). It was shown in Ref. [lj| that although e 
has non-commensurate components, the delta-functions 
merge and the conductivity may become a smooth func- 
tion of its parameters. In this paper we do not restrict 
the choice of E and e like in Ref. IKI. 



C. Effective dimension of the 0- function 

The properties of the matrix-O-function strongly de- 
pend on the number d of the nonzero components of a. 
This number can be treated as an effective dimension of 
the 8-function. Generally 1 < d < N. Without loss of 
the generality we can assume that the first d components 
of a are nonzero while the others take zero values. We 
can apply a reduction procedure if d < N: to take the 
sum over with i > d in Eq. ((HJ). We finally obtain 
a function very similar to the matrix-theta function but 
now with a summation over d-dimensional integer vec- 
tors. For example, d — 2 < N for the problem considered 
in Ref. Q and d = N for Ref. [|. 

D. Discrete environment spectra 

The problem of the evolution of a quantum system in- 
teracting with a (quantum) environment of (soft) modes 
is being studied for more than 60 years, but more im- 
portant to be solve than ever, see e.g., Refs. [TTI. Il2l. I24T- 
1301 ] - Usually it is implied that the environment has a 
continues spectrum of modes. That condition is impor- 
tant since it typically avoids the the appearance of the 
^-functions like in Eq. (JTJ) in observables of the quan- 
tum system. A discrete environment, on the other hand, 
cannot simply absorb arbitrary amounts of energy, but 
rather only discrete energy portions in quantums on the 
order of its level spacing. In other words, the set of quan- 
tum modes can serve as a "bath" only when its spectrum 
is continuous [24| , which we will clarify below. 

An important problem where a quantum system in- 
teracts with an environment, is the problem of quan- 
tum transport through an ultrasmall tunnel junction, see 
Fig. [SJ and Ref. [l2j for a review. If the contacts of the 
junction are superconducting, the supercurrent is given 
by 

I s (V)=nE 2 j(P(2V)-P(-2V)), (25) 

where -P(w) is the probability to exchange energy oj with 
the environment and Ej is the Josephson energy of the 
junction. 

Next, we concentrate on the shape of P. According to 
Ref. [ll 

P(w) = i j dt{exp[ie*4>(t)]exp[-ie*$(0)])e iut , (26) 
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a) 



R 



0^— V 




FIG. 6. A tunnel junction interacting with the environment 
of electromagnetic or many-body excitations sketched as ef- 
fective impedance Z. The wavy lines show schematically the 
exchange of the energy between the junction and the environ- 
ment [301 . b) The environment is not necessary described by 
the equilibrium density matrix [3(| • Here we sketched the en- 
vironment interacting with a thermal bath and the quantum 
system. 



where <\> is the "charge transfer" operator, e* = 2 for 
Cooper pairs and (...) denotes averaging over the den- 
sity matrix p of the environment. [We recall that we 
have chosen units where fcs = 7j=e = c=l]. Since 
the environment is isolated, it can be described by the 
Hamiltonian H cnv with a set of energies e a and eigen- 
functions \a). In general we have p = w a p\a)(/3\. Writ- 
ing the phase operator in Heiscnbcrg representation ex- 
plicitly, we get cxp[ie*0(i)] = U' 1 ^) exp[ie*<f>(0)}U(t), 
where U(t) = exp(—iH env t). This way we find: 

P(w)= 5( E a -Ep + E 7 -En + uj)x 

w Xa (e ie ^ (0) ) ^ (e- fa '* >) a . (27) 

Thus P(u)) is represented as the superposition of 8- 
functions. Only if the spectrum of the environment 
is continuous we can introduce a continuous density of 
states and integrate out the 6- functions in Eq. (f27)) . 

As long as the spectrum of the environment is discrete, 
we can assume without loss of generality that the Greek 
indices labeling the levels are a set of integer numbers. 
If the environment is represented by oscillator modes (or 
quantum rotators), E a — Ep + E n — E n + lu is linear in 
the indices labeling the environment states like in Eq. ([TJ . 
If, in addition, the environment is in local equilibrium, 
such that w\ a cx S\. a exp(—E a /T env ), where T env is the 
environment temperature, the structure of Eq. (|25p is the 
same as the structure of Eq. ((T|) . 

To complete our investigation we focus on the quasi- 
particle current J q through the tunnel junction shown in 
Fig. |U According to Refs. [U El [US] we have 



J„ 



(28) 



where F* (I 7 ) is the tunneling rate from the left (right) 
to the right (left), and, for a single junction, 



/«(l-/^)P(e-e'), 



(29) 



where / < - 1,2 - 1 are the electronic distribution functions 
within the left (right) electrodes, and R T is the bare tun- 
nel resistance, representing the interaction of electrons 
with the bath. Here P(io) is the probability for electron 
quasi-particle to lose the energy u> to the environment; it 
is given by Eq. (f27|) with e* = 1. The backward scatter- 
ing rate is given by f cx J ee , / e (2) (l - f { t P)P{e - e'). [If 
the contacts are superconducting, we can account for that 
by introducing the quasi-particle densities of the states in 
the contacts [12j in Eq. ([2"9"]l. Eq. (J3SJ would then give the 
quasi-particle current.] If the environment is absent and 
the relaxation is provided by a phonon bath, P (e) = 5(e) 
and Eq. (|2"9"|) reproduces the conventional Ohm law. It 
follows from Eqs. (f28|) - (|29 |) that the integration over the 
energy removes the (5-functions of the environment and 
the quasiparticle current is not as singular as the super- 
current when the environment has a discrete spectrum of 
modes. 

To summarize, the problem we solve in this paper is 
closely related to the problem of a discrete environment 
exchanging energy with a quantum conductor. It should 
be emphasized that the environment should not be neces- 
sary located somewhere outside the quantum conductor. 
On the contrary, it can be part of the conductor itself. 
Such a situation is realized in J J As, see e.g., Ref. or in 
very dirty conductors where the transition between dis- 
crete and continuous "built-in" environments is closely 
related to many-body localization, Refs. [30l - [32j . 



CONCLUSIONS 



We have shown that observables having a form as de- 
scribed by Eq. (JT|) can be rewritten in terms of the matrix 
theta-function. We demonstrated that the (5-functions 
typically condensate in the mesoscopic regime such that 
the observable / we focus on becomes a nearly continuous 
function. We found the envelop function for / analyti- 
cally and therefore, our results can help to understand 
and describe transport properties in a number of strongly 
correlated systems. 
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Appendix A: Parameters of the B-function. 
Analytical expressions. 

a. Identity 1 

The matrix E can be diagonalizcd by the orthogo- 
nal transformation U, such that U T \U = E, where 
A = diag(Ai, . . . , Ajv) is the diagonal matrix of the eigen- 
values. Then expanding the determinant of G over 7 we 
get the following relation: 



q = lim 7 det G = det E 

CT->0 ' ^ 



(e~,) 2 



(Al) 



e. Identity 6 

Finally, we find the inverse of the matrix G at a — > 0. 
This is somewhat tricky since G becomes singular while, 
following Eq. (|A1[) . G -1 becomes degenerate. Neverthe- 
less there is a finite nontrivial limit for G -1 at a — >• 0: 



K = lim G = cT 1 det(E) C/ T - 

CT— >0 



/ fcl 

V 



A 2 Ai 



Ajv Ai 



_eie2_ 
Ai A2 

fc 2 



ejye 2 
A/v A2 



where e — U • e. 



The diagonal elements have the str 



eiejy 
Ai Ajv 

C2GN 

A 2 A jv 



ucturc: 



17 (A7) 



b. Identities 2,3 



1 J 



(A8) 



Similarly we can prove, using an induction procedure, 
that 



where 2 = 1, 



lim 7 1 e,(G 1 ) l] e j = 1, 



(A2) 



Appendix B: Diagonal i5-matrix 



and 



p = -]im 1 - 1 [ 7 - 1 e i (G-%e j -l] 

CT— >-0 



det E 
9 

It should be noted that p > 0. 



= £ 



Ai 



We assume that e is an eigenvector of matrix E corre- 
sponding - without loss of generality - to eigenvalue Ai. 
Then e = |e|(l, 0, . . . , 0) T , and therefore 



g = V- det E 
(A3) Ai 

K = tr r diag(0,A2 1 ,A3 1 ,.. 

where Aj>i = 2tt / Vj. 

Most important in this case is 



(Bl) 

,A^)C7, (B2) 



c. Identity 4 



a = e. 



(B3) 



lim 7 1 (G 1 ) ij e : j 

CT— f0 



(A4) 



We will see that the vector a is the characteristic direc- 
tion of the i?-function. It is worth noting that generally 
e and a are not parallel in Euclidean space. The vector 
a can be found explicitly: 



( t\ 



det(.E) U T 



X a 



(A5) 



This identity can be proven in a diagonal representation 
of E. Without loss of the generality e corresponds to 
the eigenvalue Ai of E. Then we can take again e = 
|e|(l,0,...,0) T and 



C/G _1 C/ T = diag 



1 



1 



Ai +7~ 



1 

' Aat 



a= lim 1 - 1 UG- 1 U T e = e. 

<7->0 



(B4) 
(B5) 



So, a = U T a = U T e = e. 



Appendix C: numerical investigation 



d. Identity 5 



It follows from Eq. (|A4|) that 

limfG- 1 ),,^ = 0. 

CT— >0 



(A6) 



An important task is to verify the accuracy of Eq. (|2"Tj) 
and the envelope approximation of the matrix 0-function 
numerically. The required calculations of the matrix 0- 
function for N > 1 requires care since one should sum 
up many quickly oscillating functions. We have found 
an effective numerical algorithm explicitly relying on the 
existence of the soft mode e of the matrix T. 
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FIG. 7. (Color online) The sketch of the set of points m, in 
Eq. ([9|, that belong to the "cylinder" volume and satisfy the 
condition mTm < 1. Here the solid line is directed along 
e = (1.04, 0.96) r . The green diamonds distinguish the points 
of the effective crystal that give the leading contribution to 6 
for T2 < 10; the black points should be added if T2 < 2 while 
the small blue points become important for T2 < 1. 



The matrix T can be expressed explicitly through its 
eigenvalues t\ , . . . , tjv and the corresponding eigenvec- 
tors e and h^, i = 2, . . . , JV [e _L hW] 



/v 



r=T 1 e-e T + Y j nh {l) -(h«) T . 



(CI) 



i=2 



The set of m, see Eq. form a ./V-dimensional cu- 
bic crystal. We construct a cylinder in this space with 
generating lines e and an elliptical support with main di- 
rections h^K Doing numerical calculations we take into 
account only points of the crystal, m, in Eq. © that be- 
long to the cylinder volume, see Fig. [7J The proportions 
of the cylinder depend on the calculation accuracy. The 
ratio of the cylinder height to its characteristic diameter 
is of the order of min{Tj>i}/Ti S> 1. The case Tj>i ^ 1 is 
the most interesting since then the (5-functions in Eq. ([TJ) 



are expected to condensate. Then only m, the nearest 
neighbors to the line directed along e, should be taken 
into account, see green diamonds in Fig. [7] and red dots 
in Fig. f2] This property strongly reduces the numeri- 
cal efforts compared to a direct calculation of I(v) using 
Eq. Q. 

A typical graph of 0(w,T) for N = 8 is shown in 
Fig. [3] [black line], the red line shows the envelop func- 
tion represented according to Eq. (p~8|) by the one di- 
mensional (Jacoby) O-function. The inset shows that 
the ^-function width a is the minimal characteristic "fre- 
quency" of the matrix 9-function. It follows also that 
the matrix 0-function can be approximated by the one- 
dimensional (Jacoby) 0-function with specially chosen 
parameters. For producing the graph [3] we have used T 
from Eq. (|14[) . This minimal model we have chosen for 
two reasons: 1) T is parameterized by a minimal set of 
parameters, 2) this model has some relation to physical 
applications, see Refs. 0,0. 

The vector e we expressed as e = e'°) + Se, where 
we chose e^ ^ = (1, 1, . . . , 1) according to Rcf. [10| . The 
components of 5e are generated by a Gaussian random 
number generator with zero average and the variance a. 
For simplicity we used the restriction, ^ i Sei = 0. Then 

we get r* = 27T ^- + a 2 that agrees with Ref. ■ 

If e = e*- ', the components of e are commensurate 
and 0(w) behaves just as the superposition of 5 functions 
shifted by a constant period of the order unity along the 
z-axis as follows from Eq. (QJ. But if Se =^ then the 
components of e are not commensurate, there is no peri- 
odicity in the (5-function set, as it is illustrated in Fig. [3J 
for a = 0.05. 

We implied nearly everywhere above that j3 < 1. If this 
is not the case and j3 > 1, P n in Eq. ([TJ) quickly decay 
with growing |n| and we can disregard all n with |n| > 1 
within acceptable error bars. It is clear that in that case 
I(oj) [as well as $] behaves always as a set of ^-functions 
with some finite distance from each other. Disorder in e 
only slightly shifts the positions of the (5-functions on the 
w-axis and there is no (5-function condensation. What 
happens when we go from /3 < 1 to f3 > 1 is illustrated 
in Fig. 01 When f3 > 1 the theta function behaves as 
a discrete set of the (5-pcaks. With decreasing j3 more 
peaks develop, for /3wl the peaks start merging and for 
/?< 1 the matrix theta-function becomes nearly smooth. 
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